#-----------------------------------#
# Alternative community detection   #                  
#-----------------------------------#
rm(list = ls())


## Function for later
range01 <- function(x){(x-min(x, na.rm = T))/(max(x, na.rm = T)-min(x, na.rm = T))}


## data sets
adjacency_matricies <- readRDS("intl_order/adjacency_matricies_replication.rds")
war.data <- readRDS("data/war_data_replication.rds")

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- max_year$max_year[i] + 1
  temp_data <- subset(community.blondel, names %in% names_temp & year == year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(subset(temp_data, group == to[z])$names)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

overlap_coef_louvain <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")

war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(group.x != group.y, 1, 0)
)

overlap_coef_louvain_cross_sectional <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")
summary(overlap_coef_louvain_cross_sectional)

set.seed(1)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_leiden(g, weights = E(g)$weight, n_iterations = 1000, objective_function = "modularity")
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- max_year$max_year[i] + 1
  temp_data <- subset(community.blondel, names %in% names_temp & year == year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(subset(temp_data, group == to[z])$names)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_leiden(g, weights = E(g)$weight, n_iterations = 1000, objective_function = "modularity")
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

orig <- readRDS("intl_order/international_order_replication.rds")
community.leiden <- community.blondel

leiden_check <- dplyr::select(
  community.leiden, 
  c(
    names, year, interlayer
  )
) %>% dplyr::rename(
  leiden_coms = interlayer
) %>% distinct(
  names, year, .keep_all = T
)

check_data <- left_join(
  orig, 
  leiden_check
)

saveRDS(check_data, file = "intl_order/international_order_leiden.rds")

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

overlap_coef_leiden <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")

war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(group.x != group.y, 1, 0)
)

overlap_coef_leiden_cross_sectional <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")
summary(overlap_coef_louvain_cross_sectional)

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- max_year$max_year[i] + 1
  temp_data <- subset(community.blondel, names %in% names_temp & year == year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)

war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

leiden_ts <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")

war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(group.x != group.y, 1, 0)
)

leiden_cross_sectional <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")
summary(leiden_ts)
summary(leiden_cross_sectional)

community.blondel <- readRDS("intl_order/international_order_replication.rds")


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(group.x != group.y, 1, 0)
)

louvain_cross_sectional <- glm(fatal_mid ~ out_group + contig + democ_joint  + autoc_joint + cinc_ratio   + major_power_joint  + poly(year, 2), data = war.data.temp.test, family = "binomial")

summary(louvain_cross_sectional)


extract.glm <- function(model,
                        include.deviance = TRUE,
                        include.nobs = TRUE,
                        ...) {
  s <- summary(model)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names <- rownames(s$coefficients)
  coefficient.names[coefficient.names == "IIntercept)"] <- "Constant"
  coefficient.names[coefficient.names == "out_group"] <- "Out-cooperative regime relation"
  coefficient.names[coefficient.names == "democ_joint"] <- "Joint democracy"
  coefficient.names[coefficient.names == "autoc_joint"] <- "Joint autocracy"
  coefficient.names[coefficient.names == "major_power_joint"] <- "Joint major powers"
  coefficient.names[coefficient.names == "cinc_ratio"] <- "Cinc ratio"
  coefficient.names[coefficient.names == "contig"] <- "Contiguity"
  
  co <- s$coefficients[, 1]
  se <- s$coefficients[, 2]
  pval <- s$coefficients[, 4]
  
  
  gof <- numeric()
  gof.names <- character()
  gof.decimal <- logical()
  
  if (include.deviance == TRUE) {
    gof <- c(gof, deviance(model))
    gof.names <- c(gof.names, "Deviance")
    gof.decimal <- c(gof.decimal, TRUE)
  }
  if (include.nobs == TRUE) {
    n <- s$df[2]
    gof <- c(gof, n)
    gof.names <- c(gof.names, "Num. obs.")
    gof.decimal <- c(gof.decimal, FALSE)
  }
  
  
  
  tr <- createTexreg(
    coef.names = coefficient.names,
    coef = co,
    se = se,
    pvalues = pval,
    gof.names = gof.names,
    gof = gof,
    gof.decimal = gof.decimal
  )
  return(tr)
}

texreg_louvain_cross_sectional              <- extract.glm(louvain_cross_sectional)
texreg_overlap_coef_louvain                 <- extract.glm(overlap_coef_louvain)

texreg_leiden_cross_sectional               <- extract.glm(leiden_cross_sectional)
texreg_leiden_ts                            <- extract.glm(leiden_ts)
texreg_overlap_coef_leiden                  <- extract.glm(overlap_coef_leiden)
texreg_overlap_coef_leiden_cross_sectional  <- extract.glm(overlap_coef_leiden_cross_sectional)

texreg(list(texreg_louvain_cross_sectional, texreg_overlap_coef_louvain,  
            texreg_leiden_cross_sectional, texreg_overlap_coef_leiden), stars = 0.05, digits = 2)

alt_community_mods <- texreg(list(texreg_louvain_cross_sectional, texreg_overlap_coef_louvain,  
            texreg_leiden_cross_sectional, texreg_overlap_coef_leiden), stars = 0.05, digits = 2)
write.table(alt_community_mods, file = "data_visualization/si_table_10.txt")

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

two_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(two_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

three_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(three_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

four_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(four_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))


set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

five_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(five_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))


set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5, max_year$max_year[i] + 6)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

six_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(six_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5, max_year$max_year[i] + 6, max_year$max_year[i] + 7)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

seven_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(seven_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5, max_year$max_year[i] + 6, max_year$max_year[i] + 7, max_year$max_year[i] + 8)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

eight_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(eight_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5, max_year$max_year[i] + 6, max_year$max_year[i] + 7, max_year$max_year[i] + 8, max_year$max_year[i] + 9)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

nine_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(nine_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_louvain(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1, max_year$max_year[i] + 2, max_year$max_year[i] + 3, max_year$max_year[i] + 4, max_year$max_year[i] + 5, max_year$max_year[i] + 6, max_year$max_year[i] + 7, max_year$max_year[i] + 8, max_year$max_year[i] + 9, max_year$max_year[i] + 10)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_louvain(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

ten_year_lag <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(ten_year_lag, cluster = c("ccode1", "ccode2", "dyad", "year"))

diff_lags <- data.frame(
  lag = c(
    paste0("Lag of ", 2:10, " years")
  ),
  est = c(
    two_year_lag$coeftable$Estimate[1],
    three_year_lag$coeftable$Estimate[1],
    four_year_lag$coeftable$Estimate[1],
    five_year_lag$coeftable$Estimate[1],
    six_year_lag$coeftable$Estimate[1],
    seven_year_lag$coeftable$Estimate[1],
    eight_year_lag$coeftable$Estimate[1],
    nine_year_lag$coeftable$Estimate[1],
    ten_year_lag$coeftable$Estimate[1]
  ),
  se = c(
    summary(two_year_lag)$se[1],
    summary(three_year_lag)$se[1],
    summary(four_year_lag)$se[1],
    summary(five_year_lag)$se[1],
    summary(six_year_lag)$se[1],
    summary(seven_year_lag)$se[1],
    summary(eight_year_lag)$se[1],
    summary(nine_year_lag)$se[1],
    summary(ten_year_lag)$se[1]
  )
) %>% mutate(
  lb = est - 2 * se,
  ub = est + 2 * se,
  lag = factor(lag, levels = paste0("Lag of ", 2:10, " years")),
  sig = ifelse(lag == "Lag of 7 years", "Not statistically\nsignificant", "Statistically\nsignificant"),
  facet_lab = "Varying the year lag for interlayer detection"
)


diff_lag_plot <- ggplot(diff_lags, aes(x = est, xmin = lb, xmax = ub, col = sig, y = lag)) + 
  geom_errorbar(width = 0.25, size = 1.25) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.25) +
  scale_colour_brewer(palette = "Set1") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 16, color = "black"),
        strip.text.x = element_text(size = 18, color = "black"), 
        axis.text = element_text(size = 18,  color = "black"), 
        axis.title = element_text(size = 18,  color = "black"),
        axis.ticks.x = element_blank()) + 
  labs(x = "", y = "") +
  facet_wrap(~facet_lab)


ggsave("data_visualization/different_lags.png", diff_lag_plot, width = 10, height = 9, units = "in")  


set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_walktrap(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

#diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_walktrap(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

walktrap_mod <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(walktrap_mod, cluster = c("ccode1", "ccode2", "dyad", "year"))

set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_fast_greedy(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

#diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_fast_greedy(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)

fast_greedy_model <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(fast_greedy_model, cluster = c("ccode1", "ccode2", "dyad", "year"))


set.seed(9732)     
community.blondel <- lapply(adjacency_matricies, function(x){ 
  g <- graph_from_adjacency_matrix(x, mode = "undirected", weighted = T)
  comm <- cluster_label_prop(g, weights = E(g)$weight)
  glag.info <- data.frame(
    names = comm$names,
    group = comm$membership
  ) 
  return(glag.info)
}) 

for (i in 1:length(community.blondel))
{
  community.blondel[[i]]$year <- i + 1815 - 1
  if (i != 1)
  {
    community.blondel[[i]]$group <- community.blondel[[i]]$group + max(community.blondel[[i - 1]]$group) 
  }
}

community.blondel <- do.call(bind_rows, community.blondel)
community.blondel$group <- as.numeric(as.factor(community.blondel$group))
summary(community.blondel$group)  


group_matrix <- matrix(0, nrow = length(unique(community.blondel$group)), ncol = length(unique(community.blondel$group)))
colnames(group_matrix) <- rownames(group_matrix) <- unique(community.blondel$group)

gc()    
max_year <- plyr::ddply(
  community.blondel, 
  ~group, 
  summarize, 
  max_year = max(year)
) %>% arrange(
  max_year
)

group_names <- unique(colnames(group_matrix))
edgelist <- data.frame(from = NA, to = NA, weight = NA)

for (i in 1:nrow(max_year))
{
  names_temp <- subset(community.blondel, group == max_year$group[i])$names
  year_temp  <- c(max_year$max_year[i] + 1)
  temp_data <- subset(community.blondel, names %in% names_temp & year %in% year_temp)
  from <- max_year$group[i]
  to <- unique(temp_data$group)
  weights <- c()
  
  if (length(to) > 0)
  {
    ## Jaccard similarity for each layer
    for (z in 1:length(to))
    {
      union_data <- union(names_temp, subset(temp_data, group == to[z])$names)
      sender_vec <- ifelse(union_data %in% names_temp, 1, 0)
      reciever_vec <- ifelse(union_data %in% subset(temp_data, group == to[z])$names, 1, 0)
      weights[z] <- length(intersect(names_temp, subset(temp_data, group == to[z])$names))/length(union_data)
    }
    
    edgelist <- bind_rows(
      edgelist,
      data.frame(
        from = rep(from, length(to)),
        to = to,
        weight = weights
      )
    ) %>% filter(
      !is.na(from)
    )
  }
}
edgelist$weight <- round(edgelist$weight, 2)    
for (i in 1:nrow(edgelist))
{
  group_matrix[rownames(group_matrix) == edgelist$from[i],
               colnames(group_matrix) == edgelist$to[i]] <- edgelist$weight[i] 
  
  group_matrix[rownames(group_matrix) == edgelist$to[i], 
               colnames(group_matrix) == edgelist$from[i]] <- edgelist$weight[i] 
}

#diag(group_matrix) <- 2
group_matrix[group_matrix < 0] <- 0

gc()    
g <- graph_from_adjacency_matrix(group_matrix, mode = "undirected", weighted = T)
set.seed(477)
comm <- cluster_label_prop(g, weights = E(g)$weight)
between.group <- data.frame(
  group = comm$names,
  interlayer = comm$membership
)



community.blondel <- community.blondel[,!(colnames(community.blondel) %in% c("interlayer", "n"))]
community.blondel <- mutate(
  community.blondel,
  group = as.character(group)
) %>% left_join(
  between.group,
  by = c(
    "group"
  )
) 

group.size <- mutate(
  community.blondel,
  n = 1
) %>% aggregate(
  n ~ year + interlayer, data = ., FUN = sum
)

community.blondel <- left_join(
  community.blondel,
  group.size,
  by = c(
    "year", "interlayer"
  )
)


war.data.temp.test <- distinct(
  war.data, 
  namea, nameb, year, mid, fatal_mid, .keep_all = T
) %>% left_join(
  community.blondel,
  by = c(
    "namea" = "names", "year"
  )
) %>% left_join(
  community.blondel,
  by = c(
    "nameb" = "names", "year"
  )
) %>% mutate(
  dyad = paste(namea, nameb, sep = "-"),
  contig = if_else(contig %in% 1:2, 1, 0),
  out_group = if_else(interlayer.x != interlayer.y, 1, 0)
)


cluster_label_prop_model <- fixest::feglm(fatal_mid ~ out_group | ccode1 + ccode2 + dyad + year, data = war.data.temp.test, family = "binomial")
summary(cluster_label_prop_model, cluster = c("ccode1", "ccode2", "dyad", "year"))


diff_com <- data.frame(
  lag = c(
    "Walktrap", "Fast\n& greedy", "Label\npropogation"
  ),
  est = c(
    walktrap_mod$coeftable$Estimate[1],
    fast_greedy_model$coeftable$Estimate[1],
    cluster_label_prop_model$coeftable$Estimate[1]
  ),
  se = c(
    summary(walktrap_mod)$se[1],
    summary(fast_greedy_model)$se[1],
    summary(cluster_label_prop_model)$se[1]
  )
) %>% mutate(
  lb = est - 2 * se,
  ub = est + 2 * se,
  sig = "Statistically\nsignificant",
  facet_lab = "Alterantive community detection algorithms"
)

diff_comm_plot <- ggplot(diff_com, aes(x = est, xmin = lb, xmax = ub, col = sig, y = lag)) + 
  geom_errorbar(width = 0.10, size = 1.25) + 
  theme_bw() + 
  geom_vline(xintercept = 0, linetype = "dashed", size = 1.25) +
  scale_color_manual(values = "#377EB8") + 
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        legend.text = element_text(size = 16, color = "black"),
        strip.text.x = element_text(size = 18, color = "black"), 
        axis.text = element_text(size = 18,  color = "black"), 
        axis.title = element_text(size = 18,  color = "black"),
        axis.ticks.x = element_blank()) + 
  labs(x = "", y = "") +
  facet_wrap(~facet_lab) + 
  scale_x_continuous(limits = c(-0.1, 1.3))


both_specifications <- cowplot::plot_grid(diff_comm_plot, diff_lag_plot)
ggsave("data_visualization/different_lags_comms.png", both_specifications, width = 17, height = 9, units = "in")  

